Estadísticas espacio-temporales y Modelado de la distribución biogeográfica de Jaguares y Pumas en Centro y Sudamerica
Autor: Jose Francisco Nuñez Obando | Posgrado Analista de Datos
26 January 2020
1 Introducción
El siguiente documento corresponde a un trabajo de investigación del curso Visualización de datos y análisis estadístico del Posgrado Analista de Datos de la Universidad Cenfotec.
2 Tema de investigación
“Estadisticas y Modelado de la distribución biogeográfica de Pumas y Jaguares de un sector del continente americano”
3 Fuentes de datos
https://www.inaturalist.org/
https://biodiversityinformatics.amnh.org/
4 Objetivos
-
- Estadísticas espacio-temporales sobre los pumas y jaguares de un sector del continente americano, en las que se evalua: frecuencia de observaciones por mes, patrones de distribución dentro del area de estudio asi como densidad por especie.
- 2.Analisis de correlacion con las variables ambientales en relacion la ocurrencia de las especies analizadas.
-
- Modelado de la distribución de los pumas y jaguares del área de estudio, segun las variables ambientales en la que estas especies se encuentran.
5 Librerias
library(dplyr)
library(tibble)
library(tidyr)
library(sp)
library(DT)
library(leaflet)
library(leaflet.extras)
library(raster)
library(rasterVis)
library(rgdal)
library(ggplot2)
library(plotly)
library(LaCroixColoR)
library(sf)
library(here)
library(devtools)
library(proto)
library(highcharter)
library(rgeos)
library(velox)
library(usdm)
library(gtools)
library(corrplot)
library(Hmisc)
library(dismo)6 Procesos de extracción, transformación y carga de las Bases de datos
En esta sección se realizan una serie de procesos que conllevan al desarrollo de las bases de datos finales
6.1 Base de datos sobre felinos
Carga de la base de datos y de las variables que la conforman
df <- read.csv("E:/Analista de datos/DAT004/Proyecto Final/Datasets/samples/wildcats_df.csv", stringsAsFactors=F, sep=",", na.strings=c("","NA"), encoding="UTF-8")
colnames(df)## [1] "X.U.FEFF.common_name" "scientific_name" "latitude"
## [4] "longitude" "created_at" "updated_at"
## [7] "place_country_name"Cambio de nombre de las columnas
df <- df[,c(1:5,7)]
names(df)[1] <- "Nombre_Comun"
names(df)[2] <- "Nombre_Cientifico"
names(df)[3] <- "Latitud"
names(df)[4] <- "Longitud"
names(df)[5] <- "Fecha_Registro"
names(df)[6] <- "Pais"
colnames(df)## [1] "Nombre_Comun" "Nombre_Cientifico" "Latitud"
## [4] "Longitud" "Fecha_Registro" "Pais"6.1.1 Cantidad de registros en la base de datos
Se imprime un resumen de datos por especie de modo que se pueda obtener un conteo inicial
df%>%
dplyr::group_by(Nombre_Comun)%>%
dplyr::summarise(length(Nombre_Comun))Se realiza una correccion de datos en la columna de “Nombre_Comun” de forma que se estandaricen algunos registros que figuran como duplicados
ndf <- df%>%
mutate(Nombre_Comun=ifelse(Nombre_Comun=="Puma de América del Norte", "Puma",Nombre_Comun))%>%
mutate(Nombre_Comun=ifelse(Nombre_Comun=="Puma de América del Sur", "Puma",Nombre_Comun))Nuevamente se imprime un resumen final de datos por especie
count_sp <- ndf%>%
group_by(Nombre_Comun)%>%
summarise(Count=length(Nombre_Comun))
plot_ly(data=count_sp, x = ~Nombre_Comun, y = ~Count, type = 'bar', text = ~Count, textposition = 'auto',
marker = list(color = c('rgba(204,204,204,1)', 'rgba(222,45,38,0.8)'),
line = list(color = c('rgba(204,204,204,1)', 'rgba(222,45,38,0.8)'),
width = 1.5))) %>%
layout(title = "Cantidad de especies registrados en el dataframe inicial",
xaxis = list(title = "Especies"),
yaxis = list(title = "Cantidad"))6.1.2 Mapa de la distribicion de las especies a nivel mundial
cf <- colorFactor(c("#ffa500", "#13ED3F"), domain=ndf$Nombre_Comun)
m <- leaflet(df) %>% addTiles(group = "OSM (default)") %>%
addProviderTiles(providers$Stamen.Toner, group = "Toner") %>%
addProviderTiles(providers$Stamen.TonerLite, group = "Toner Lite") %>%
setView(-63.54105072, -8.88123188, zoom = 2) %>%
addCircles(~Longitud, ~Latitud, popup=ndf$Nombre_Comun, group = ndf$Nombre_Comun ,weight = 3, radius=4,
color=cf(ndf$Nombre_Comun), stroke = TRUE, fillOpacity = 0.5) %>%
addLegend("bottomright", colors= c("#ffa500", "#13ED3F"), labels=c("Puma", "Jaguar"), title="Especies") %>%
addLayersControl(
baseGroups = c("OSM (default)", "Toner", "Toner Lite"),
overlayGroups = ndf$Nombre_Comun,
options = layersControlOptions(collapsed = TRUE))%>%
suspendScroll()
m6.2 Bases de datos sobre variables ambientales del área de estudio
El area de estudio con la que se pretende trabajar corresponde a Centro y Sudamrica, asi como el Caribe.
Algunas de las coberturas que se utilizan para conocer el habitat de las especies son variables climáticas que se derivan de los datos proporcionados por el Panel Intergubernamental sobre Cambio Climático y se produjeron utilizando Interpolación de basada en lecturas tomadas en estaciones meteorológicas de todo el mundo desde 1961 hasta 1990.
Los coberturas corresponden a:
*cld6190_ann : cobertura de nubes, anual
*dtr6190_ann : rango de temperatura diurna, anual
*ecoreg : Ecoregiones
*frs6190_ann : frecuencia de heladas, anual
*h_dem : Modelo Digital de Elevacion
*pre6190_ann : precipitacion, anual
*pre6190_I1 : precipitacion, enero
*pre6190_I4 : precipitacion, abril
*pre6190_I7 : precipitacion, julio
*pre6190_I10 : precipitacion, octubre
*tmn6190_ann : temperatura promedio, anual
*tmp6190_ann : temperatura minima, anual
*tmx6190_ann : temperatura maxima, anual
*vap6190_ann : presion de vapor, anual
6.2.1 Datos Ambientales
r_files <- list.files(here::here("Datasets", "coverages"),
full.names = T)
r_files## [1] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/cld6190_ann.asc"
## [2] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/dtr6190_ann.asc"
## [3] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/ecoreg.asc"
## [4] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/frs6190_ann.asc"
## [5] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/h_dem.asc"
## [6] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/pre6190_ann.asc"
## [7] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/pre6190_l1.asc"
## [8] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/pre6190_l10.asc"
## [9] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/pre6190_l4.asc"
## [10] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/pre6190_l7.asc"
## [11] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/tmn6190_ann.asc"
## [12] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/tmp6190_ann.asc"
## [13] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/tmx6190_ann.asc"
## [14] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/vap6190_ann.asc"6.2.1.1 Ploteo de las variables ambientales
st <- raster::stack(r_files)
raster::plot(st)
6.2.1.2 Muestra de una de las variables ambientales: Ecoregiones
1 Mediterranean Forests, Woodlands & Scrub 2 Tropical & Subtropical Dry Broadleaf Forests 3 Temperate Broadleaf & Mixed Forests 4 Tropical & Subtropical Coniferous Forests 5 Temperate Grasslands, Savannas & Shrublands 6 Mangroves 7 ( ) 8 Montane Grasslands & Shrublands 9 Tropical & Subtropical Grasslands, Savannas & Shrublands 10 Tropical & Subtropical Moist Broadleaf Forests 11 Flooded Grasslands & Savannas 12 Montane Grasslands & Shrublands 13 Magellanic subpolar forests 14 Rock and Ice 15 ( )
ecoreg = raster(r_files[3])
raster::plot(ecoreg, main="Ecoregiones del area de estudio")
raster::hist(ecoreg, main="Histograma de las Ecoregiones del area de estudio", col="green")
6.3 Delimitacion del area de estudio
A continuacion, se carga el un archivo espacial de tipo vectorial correspondiente a los paises que conforman el area de estudio.
area_estudio <- sf::st_read(dsn="E:/Analista de datos/DAT004/Proyecto Final/Datasets/samples/map.geojson", layer="map", quiet=TRUE)6.3.0.1 Mapa de paises del area de estudio
map.area <- plot_geo(area_estudio, name=~pais, color = I("#89C5DA"))%>% hide_legend()
map.area6.3.0.2 Mapeo de los ubicaciones registradas de las especies en el área de estudio
De previo se realiza una selección de la base de datos de las especies respecto del área de estudio
coords <- data.frame(x=ndf$Longitud, y=ndf$Latitud)
prj <- CRS("+proj=longlat +datum=WGS84 +no_defs")
points <- SpatialPointsDataFrame(coords, ndf, proj4string = prj)
wildcats <- crop(points,ecoreg)
wildcats <- data.frame(wildcats)
wildcats <- wildcats[,1:6]Para luego visualizar los puntos de observacion de las especies a estudiar
cf2 <- colorFactor(c("#ffa500", "#13ED3F"), domain=wildcats$Nombre_Comun)
m2 <- leaflet(wildcats) %>% addTiles(group = "OSM (default)") %>%
addProviderTiles(providers$Stamen.Toner, group = "Toner") %>%
addProviderTiles(providers$Stamen.TonerLite, group = "Toner Lite") %>%
setView(-63.54105072, -8.88123188, zoom = 2) %>%
addCircles(~Longitud, ~Latitud, popup=wildcats$Nombre_Comun, group = wildcats$Nombre_Comun ,weight = 3, radius=4,
color=cf2(wildcats$Nombre_Comun), stroke = TRUE, fillOpacity = 0.5) %>%
addLegend("bottomright", colors= c("#ffa500", "#13ED3F"), labels=c("Puma", "Jaguar"), title="Especies") %>%
addLayersControl(
baseGroups = c("OSM (default)", "Toner", "Toner Lite"),
overlayGroups = wildcats$Nombre_Comun,
options = layersControlOptions(collapsed = TRUE))%>%
suspendScroll()
m26.3.0.3 Visualizacion de la base de datos
Primero se realiza una selección final de la base de datos e inclusion de nuevas variables
wildcats$Hora_Registro <- strftime(wildcats$Fecha_Registro, format="%H") #col7
wildcats$Fecha_Registro <- as.POSIXct(wildcats$Fecha_Registro, tz="", format="%Y-%m-%d") #col5
wildcats$Y.M.R <- strftime(wildcats$Fecha_Registro, format="%Y-%m") #col8
wildcats$M.R <- format(wildcats$Fecha_Registro, "%m") #col9
wildcats$Y.R <- format(wildcats$Fecha_Registro, "%Y") #col10
wildcats <- wildcats[, colnames(wildcats)[c(1:4,6,5,8,10,9,7)]] #Reorganizacion de las columnasSe procede a visualizar la base de datos
datatable(wildcats, extensions = c('Buttons','FixedColumns', 'RowReorder'),
filter = list(position = 'top', clear = FALSE),
options = list(dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print'), dom = 't',
scrollX = TRUE,
fixedColumns = TRUE,
rowReorder = TRUE, order = list(c(0 , 'asc'))))7 Resultados
7.1 Estadísticas espacio-temporales de las especies estudiadas
En este primer apartado se realizan analisis sobre estadísticas espacio-temporales respecto de los pumas y jaguares de un sector del continente americano
7.1.1 Cantidad de especies segun ubicacion y temporalidad
7.1.1.1 Especies por pais y mes de observacion
plot_ly(x= wildcats$Pais,split=wildcats$Nombre_Comun, color= I(cf2(wildcats$Nombre_Comun)),
frame=wildcats$M.R, type = "histogram", alpha = 0.6)%>%
layout(barmode = "overlay")7.1.1.2 Cantidad de Jaguares por pais
c.jaguars <- wildcats%>%
dplyr::filter(Nombre_Comun == "Jaguar")%>%
dplyr::group_by(Pais)%>%
dplyr::summarise(Cant=n())
c.pumas <- wildcats%>%
dplyr::filter(Nombre_Comun == "Puma")%>%
dplyr::group_by(Pais)%>%
dplyr::summarise(Cant=n())data(worldgeojson, package = "highcharter")
highchart() %>%
hc_title(text = "Cantidad de Jaguares") %>%
hc_subtitle(text = "por pais") %>%
hc_add_series_map(map = worldgeojson, df = c.jaguars, name = "cantidad de individuos",
value = "Cant", joinBy = c("name", "Pais"), allAreas=FALSE, animation=TRUE)%>%
hc_colorAxis(min= 1,
type='logarithmic',
stops = color_stops())%>%
hc_tooltip(useHTML = TRUE, headerFormat = "",
pointFormat = "Pais: {point.name}/ Cantidad de especies: {point.Cant}")%>%
hc_legend(layout= 'horizontal',
verticalAlign= 'top')%>%
hc_mapNavigation(enabled = TRUE) 7.1.1.3 Cantidad de Pumas por pais
data(worldgeojson, package = "highcharter")
highchart() %>%
hc_title(text = "Cantidad de Pumas") %>%
hc_subtitle(text = "por pais") %>%
hc_add_series_map(map = worldgeojson, df = c.pumas, name = "cantidad de individuos",
value = "Cant", joinBy = c("name", "Pais"), allAreas= FALSE, animation=TRUE)%>%
hc_colorAxis(min= 1,
type='logarithmic',
stops = color_stops())%>%
hc_tooltip(useHTML = TRUE, headerFormat = "",
pointFormat = "Pais: {point.name}/ Cantidad de especies: {point.Cant}")%>%
hc_legend(layout= 'horizontal',
verticalAlign= 'top')%>%
hc_mapNavigation(enabled = TRUE) 7.1.2 Cantidad de individuos observados en el tiempo
7.1.2.1 Especies por mes y año
c.jaguars <- wildcats%>%
dplyr::filter(Nombre_Comun == "Jaguar")%>%
dplyr::group_by(Fecha_Registro)%>%
dplyr::summarise(Cant=n())
c.pumas <- wildcats%>%
dplyr::filter(Nombre_Comun == "Puma")%>%
dplyr::group_by(Fecha_Registro)%>%
dplyr::summarise(Cant=n())
ct.jaguars <- ts(data=c.jaguars, start = c(2011,10), end = c(2019,12), frequency = 12)
ct.pumas <- ts(data=c.pumas, start = c(2011,10), end = c(2019,12), frequency = 12)highchart(type = "stock") %>%
hc_title(text = "Cantidad de especies") %>%
hc_subtitle(text = "en la linea de tiempo de observacion") %>%
hc_add_series(ct.jaguars[,2], type = "line", pointInterval= 30*24*3600*1000, name="Cant de jaguares observados a la fecha", color = "#ffa500")%>%
hc_add_series(ct.pumas[,2], type = "line", pointInterval= 30*24*3600*1000, name="Cant de pumas observados a a fecha", color = "#13ED3F")%>%
hc_xAxis(title= list(text='Fechas'), type='datetime')%>%
hc_yAxis(title= list(text='Cantidad de especies observadas'), opposite = FALSE)7.1.2.2 Descomposicion y comportamiento de observaciones de los individuos en el tiempo
decomp_jaguars <- decompose(ct.jaguars[,2])
decomp_pumas <- decompose(ct.pumas[,2])
plot(decomp_jaguars)
par(adj = 1)
title(sub = "Jaguares", cex.sub = 1.5, font.sub = 3, col.sub = "#ffa500")
plot(decomp_pumas, sub="Pumas")
par(adj = 1)
title(sub = "Pumas", cex.sub = 1.5, font.sub = 3, col.sub = "#13ED3F")
7.1.3 Patron de distribucion de las especies
7.1.3.1 Distribucion espacial de las especies por mes
e.m <- ggplot(wildcats, aes(x=Longitud, y=Latitud, frame=M.R, text = paste("Pais:", Pais))) +
borders("world", xlim = c(-120, -30), ylim = c(-50, 25))+
geom_point(data=filter(wildcats), shape = 21, colour = cf2(wildcats$Nombre_Comun), size = 1)+
facet_grid(cols = vars(Nombre_Comun)) +
labs(title = "Distribucion de las especies observadas segun mes")
ggplotly(e.m)%>%
animation_opts(
2000, easing = "elastic", redraw = FALSE
)%>%
animation_slider(
currentvalue = list(prefix = "Mes ", font = list(color="red"))
)7.1.3.2 Patron geometrico de distribucion
pd <- qplot(data = wildcats, x=Longitud, y=Latitud) +
stat_ellipse(geom = "polygon", alpha = 1/2, aes(fill = Nombre_Comun, frame=M.R))+
borders("world", xlim = c(-120, -30), ylim = c(-50, 25))+
geom_point(data=filter(wildcats), aes(frame=M.R), shape = 21, fill = cf2(wildcats$Nombre_Comun), size = 2)+
facet_grid(cols = vars(Nombre_Comun)) +
labs(title = "Patron de Distribucion de las Especies")
ggplotly(pd)%>%
animation_opts(
2000, easing = "elastic", redraw = FALSE
)%>%
animation_slider(
currentvalue = list(prefix = "Mes ", font = list(color="red"))
)7.1.3.3 Densidad de cada una de las especies
dn <- ggplot(wildcats, aes(x=Longitud, y=Latitud)) +
borders("world", regions=".", xlim = c(-120, -30), ylim = c(-50, 25))+
stat_density2d(aes(fill = stat(level), frame=wildcats$M.R), geom="polygon", alpha=0.3) +
geom_point(data=filter(wildcats), aes(frame=M.R), shape = 21, fill = cf2(wildcats$Nombre_Comun), size = 2, alpha=0.3)+
facet_grid(cols = vars(Nombre_Comun)) +
scale_fill_viridis_c(option = "viridis") +
theme(legend.position = "magma") +
labs(title = "Densidad segun distribucion por especie")
ggplotly(dn)%>%
animation_opts(
4000, easing = "elastic", redraw = FALSE
)%>%
animation_slider(
currentvalue = list(prefix = "Mes ", font = list(color="red"))
)7.1.3.4 Especies segun ecoregion en la que se encuentran
Ecoregiones/Biomas 1 Mediterranean Forests, Woodlands & Scrub 2 Tropical & Subtropical Dry Broadleaf Forests 3 Temperate Broadleaf & Mixed Forests 4 Tropical & Subtropical Coniferous Forests 5 Temperate Grasslands, Savannas & Shrublands 6 Mangroves 7 ( ) 8 Montane Grasslands & Shrublands 9 Tropical & Subtropical Grasslands, Savannas & Shrublands 10 Tropical & Subtropical Moist Broadleaf Forests 11 Flooded Grasslands & Savannas 12 Montane Grasslands & Shrublands 13 Magellanic subpolar forests 14 Rock and Ice 15 ( )
ecoreg_spec <- raster::extract(ecoreg, wildcats[,4:3]) %>%
cbind(wildcats, .) %>%
as.data.frame()
ecoreg_spec <- na.omit(ecoreg_spec)
ecoreg_spec$Ecoreg <- as.character(ecoreg_spec$.)
ggplot(ecoreg_spec, aes(x=Longitud, y=Latitud, color = ecoreg_spec$Ecoreg)) +
borders("world", xlim = c(-120, -30), ylim = c(-50, 25))+
geom_point(data=filter(ecoreg_spec), shape = 21, size = 1)+
scale_colour_manual(values = lacroix_palette("PassionFruit", n=14))+
facet_grid(cols = vars(Nombre_Comun)) +
labs(title = "Especies observadas segun ecoregion en la que se encuentran")
ggplot(ecoreg_spec, aes(y=Ecoreg, x=Nombre_Comun, colour = Ecoreg)) +
geom_count(alpha=0.5) +
scale_colour_manual(values = lacroix_palette("PassionFruit", n=14))+
labs(title = "Cantidad de especies observadas por ecoregion",
x = "Especie",
y = "Ecoregiones",
size = "")
7.2 Analisis de correlacion con las variables ambientales
Proceso de union de los puntos de observacion de las especies y los raster de las variables ambientales
#Se extraen temas relevantes de las variables ambientales
select_r = dropLayer(st,3)
jp_swd <- raster::extract(select_r, wildcats[,4:3]) %>%
cbind(wildcats, .) %>%
as.data.frame()
#Se elimina las filas con valores nulos
jp_swd <- na.omit(jp_swd) 7.2.1 Ploteo de la correlacion de variables
mt <- cor(jp_swd[,11:23])
corrplot(mt, method = 'circle', type = 'upper')
7.2.2 Variables que menos se relacionan entre si
Esto se hace a traves del analisis en regresion lineal VIF, el cual es el factor de inflación de la varianza. Este cuantifica la intensidad de la multicolinealidad en un análisis de regresión normal de mínimos cuadrados. Proporciona un índice que mide hasta qué punto la varianza (el cuadrado de la desviación estándar estimada) de un coeficiente de regresión estimado se incrementa a causa de la colinealidad.
# Analisis VIF
vif.res <- vif(x = jp_swd[,11:ncol(jp_swd)])
vif.step <- vifstep(x = jp_swd[,11:ncol(jp_swd)], th = 10)#Los factores VIF mayores que 10 implican problemas graves de multicolinealidad
vrs <- vif.step@results$Variables %>% as.character()
vif.step## 4 variables from the 13 input variables have collinearity problem:
##
## pre6190_ann tmp6190_ann tmn6190_ann tmx6190_ann
##
## After excluding the collinear variables, the linear correlation coefficients ranges between:
## min correlation ( pre6190_l10 ~ pre6190_l1 ): -0.009063609
## max correlation ( vap6190_ann ~ frs6190_ann ): -0.8508242
##
## ---------- VIFs of the remained variables --------
## Variables VIF
## 1 cld6190_ann 3.025670
## 2 dtr6190_ann 1.678499
## 3 frs6190_ann 4.102199
## 4 h_dem 2.447628
## 5 pre6190_l1 1.755003
## 6 pre6190_l10 5.388244
## 7 pre6190_l4 4.074227
## 8 pre6190_l7 2.896292
## 9 vap6190_ann 5.176151vrs## [1] "cld6190_ann" "dtr6190_ann" "frs6190_ann" "h_dem" "pre6190_l1"
## [6] "pre6190_l10" "pre6190_l4" "pre6190_l7" "vap6190_ann"7.2.3 Base de datos y variables con multicolinealidad
jp_fnl <- jp_swd %>%
as_tibble %>%
dplyr::select(Nombre_Cientifico:Longitud, vrs)
head(jp_fnl)7.2.4 Correlacion entre variables ambientales con multicolinealidad
Con respecto al resultado de la correlacion se grafican las siguientes variables:
*cld6190_ann : cobertura de nubes, anual
*dtr6190_ann : rango de temperatura diurna, anual
*frs6190_ann : frecuencia de heladas, anual
*h_dem : Modelo Digital de Elevacion
*pre6190_I1 : precipitacion, enero
*pre6190_I4 : precipitacion, abril
*pre6190_I7 : precipitacion, julio
*pre6190_I10 : precipitacion, octubre
*vap6190_ann : presion de vapor, anual)
r_corr= st[[c(1,2,4,5,7,8,9,10,14)]]
raster::pairs(r_corr)
7.3 Modelado de distribucion de especies
El primer paso es extraer los valores de los predictores en las ubicaciones de los puntos
pvals <- raster::extract(st, wildcats[,4:3]) %>%
cbind(wildcats, .) %>%
as.data.frame()
presvals <- pvals[,11:24]
predictors <- st
backgr <- randomPoints(predictors, 500)
absvals <- raster::extract(predictors, backgr)
pb <- c(rep(1, nrow(presvals)), rep(0, nrow(absvals)))
sdmdata <- data.frame(cbind(pb, rbind(presvals, absvals)))
sdmdata[,'ecoreg'] = as.factor(sdmdata[,'ecoreg'])
head(sdmdata)tail(sdmdata)summary(sdmdata)## pb cld6190_ann dtr6190_ann ecoreg frs6190_ann
## Min. :0.0000 Min. :32.00 Min. : 66 10 :688 Min. : 0.00
## 1st Qu.:0.0000 1st Qu.:58.00 1st Qu.: 99 11 :259 1st Qu.: 0.00
## Median :1.0000 Median :64.00 Median :115 9 :160 Median : 1.00
## Mean :0.6736 Mean :63.82 Mean :112 8 :123 Mean : 11.65
## 3rd Qu.:1.0000 3rd Qu.:68.00 3rd Qu.:123 5 : 74 3rd Qu.: 2.00
## Max. :1.0000 Max. :83.00 Max. :175 (Other):166 Max. :214.00
## NA's :9 NA's :9 NA's : 62 NA's :9
## h_dem pre6190_ann pre6190_l1 pre6190_l10
## Min. : 0.0 Min. : 0.00 Min. : 0.00 Min. : 0.00
## 1st Qu.: 107.0 1st Qu.: 35.00 1st Qu.: 22.00 1st Qu.: 27.50
## Median : 193.0 Median : 40.00 Median : 53.00 Median : 43.00
## Mean : 490.4 Mean : 48.76 Mean : 52.81 Mean : 54.32
## 3rd Qu.: 526.0 3rd Qu.: 66.00 3rd Qu.: 75.00 3rd Qu.: 82.00
## Max. :4920.0 Max. :138.00 Max. :173.00 Max. :210.00
## NA's :57 NA's :9 NA's :9 NA's :9
## pre6190_l4 pre6190_l7 tmn6190_ann tmp6190_ann
## Min. : 0.00 Min. : 0.00 Min. :-104.0 Min. : 35.0
## 1st Qu.: 22.00 1st Qu.: 6.00 1st Qu.: 140.0 1st Qu.:226.0
## Median : 34.00 Median : 28.00 Median : 165.0 Median :254.0
## Mean : 45.28 Mean : 42.45 Mean : 149.1 Mean :233.1
## 3rd Qu.: 65.00 3rd Qu.: 71.00 3rd Qu.: 187.5 3rd Qu.:264.0
## Max. :148.00 Max. :194.00 Max. : 227.0 Max. :282.0
## NA's :9 NA's :9 NA's :9 NA's :9
## tmx6190_ann vap6190_ann
## Min. :122.0 Min. : 10.0
## 1st Qu.:308.0 1st Qu.:204.5
## Median :324.0 Median :257.0
## Mean :312.1 Mean :229.1
## 3rd Qu.:335.0 3rd Qu.:267.0
## Max. :359.0 Max. :308.0
## NA's :9 NA's :97.3.1 Ajuste del modelo
Para ello se utiliza la función (glm) la cual permite ajustar modelos lineales generalizados. Esta funcion devuelve un objeto modelo.
En un primer modelo se trabajara con:
*cld6190_ann : cobertura de nubes, anual
*dtr6190_ann : rango de temperatura diurna, anual
*frs6190_ann : frecuencia de heladas, anual
*h_dem : Modelo Digital de Elevacion
*pre6190_I1 : precipitacion, enero
*pre6190_I4 : precipitacion, abril
*pre6190_I7 : precipitacion, julio
*pre6190_I10 : precipitacion, octubre
*vap6190_ann : presion de vapor, anual)
En el segundo con todas las variables
m1 <- glm ( pb ~ cld6190_ann + dtr6190_ann + frs6190_ann + h_dem + pre6190_l1 + pre6190_l10 + pre6190_l4 + pre6190_l7 + vap6190_ann, data = sdmdata )
class ( m1 )## [1] "glm" "lm"summary ( m1 )##
## Call:
## glm(formula = pb ~ cld6190_ann + dtr6190_ann + frs6190_ann +
## h_dem + pre6190_l1 + pre6190_l10 + pre6190_l4 + pre6190_l7 +
## vap6190_ann, data = sdmdata)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0999 -0.4707 0.1940 0.3046 1.0077
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.382e-01 1.602e-01 -3.985 7.08e-05 ***
## cld6190_ann 7.567e-03 1.927e-03 3.927 9.00e-05 ***
## dtr6190_ann 4.495e-03 7.680e-04 5.853 5.94e-09 ***
## frs6190_ann 2.029e-03 6.587e-04 3.080 0.002107 **
## h_dem -8.450e-05 2.258e-05 -3.743 0.000189 ***
## pre6190_l1 1.509e-03 4.505e-04 3.348 0.000834 ***
## pre6190_l10 4.674e-03 5.162e-04 9.054 < 2e-16 ***
## pre6190_l4 -5.923e-03 5.860e-04 -10.108 < 2e-16 ***
## pre6190_l7 1.351e-03 4.270e-04 3.164 0.001591 **
## vap6190_ann 9.905e-04 3.944e-04 2.512 0.012126 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.1829705)
##
## Null deviance: 320.40 on 1472 degrees of freedom
## Residual deviance: 267.69 on 1463 degrees of freedom
## (59 observations deleted due to missingness)
## AIC: 1690.4
##
## Number of Fisher Scoring iterations: 2m2 = glm ( pb ~ . , data = sdmdata )
m2##
## Call: glm(formula = pb ~ ., data = sdmdata)
##
## Coefficients:
## (Intercept) cld6190_ann dtr6190_ann ecoreg2 ecoreg3 ecoreg4
## 0.8481994 0.0003222 0.0076312 -0.2080109 0.2056482 -0.0591186
## ecoreg5 ecoreg6 ecoreg8 ecoreg9 ecoreg10 ecoreg11
## 0.1378429 -0.0442762 0.0486044 -0.0145809 0.0721378 0.5030165
## ecoreg12 ecoreg13 ecoreg14 frs6190_ann h_dem pre6190_ann
## 0.2519882 0.3957019 0.7397254 -0.0017491 -0.0001387 0.0124540
## pre6190_l1 pre6190_l10 pre6190_l4 pre6190_l7 tmn6190_ann tmp6190_ann
## -0.0032611 0.0012264 -0.0065611 -0.0013567 -0.0004383 0.0071411
## tmx6190_ann vap6190_ann
## -0.0082765 -0.0011003
##
## Degrees of Freedom: 1467 Total (i.e. Null); 1442 Residual
## (64 observations deleted due to missingness)
## Null Deviance: 319.2
## Residual Deviance: 216 AIC: 14077.3.2 Modelado de envoltura climática
El algoritmo BIOCLIM calcula la similitud de una ubicación comparando los valores de las variables ambientales en cualquier ubicación con una distribución porcentual de los valores en ubicaciones conocidas de ocurrencia. Cuanto más cerca del percentil 50 (la mediana), más adecuada es la ubicación. Las colas de la distribución no se distinguen, es decir, el percentil 10 se trata como equivalente al percentil 90.
Para este caso solo se usa datos de presencia, por lo que usamos ‘presvals’ en lugar de ‘sdmdata’.
bc <- bioclim ( presvals [, c ( 'cld6190_ann', 'dtr6190_ann', 'frs6190_ann', 'h_dem', 'pre6190_l1', 'pre6190_l10', 'pre6190_l4', 'pre6190_l7', 'vap6190_ann' )])
class ( bc )## [1] "Bioclim"
## attr(,"package")
## [1] "dismo"bc## class : Bioclim
##
## variables: cld6190_ann dtr6190_ann frs6190_ann h_dem pre6190_l1 pre6190_l10 pre6190_l4 pre6190_l7 vap6190_ann
##
##
## presence points: 1002
## cld6190_ann dtr6190_ann frs6190_ann h_dem pre6190_l1 pre6190_l10 pre6190_l4
## 1 66 118 0 107 72 30 34
## 2 57 99 0 31 113 110 88
## 3 60 121 1 250 14 45 14
## 5 63 123 10 757 86 43 32
## 7 56 84 1 51 48 98 25
## 9 60 119 1 275 15 46 14
## 10 58 97 0 2 121 111 87
## 12 60 95 0 2 23 58 12
## 14 66 118 0 110 72 30 34
## 15 53 116 1 200 47 130 43
## pre6190_l7 vap6190_ann
## 1 6 260
## 2 150 269
## 3 42 266
## 5 6 214
## 7 134 268
## 9 43 266
## 10 165 278
## 12 47 281
## 14 6 260
## 15 171 238
## (... ... ...)pairs ( bc )
7.3.3 Predicción del modelo
Todos estos objetos “modelo”, independientemente de su clase exacta, pueden utilizarse con la función de predicción para hacer predicciones para cualquier combinación de valores de las variables independientes. Hacer tales predicciones para algunos entornos puede ser muy útil para explorar y comprender las predicciones del modelo. A continuación, se realizan predicciones con el objeto del modelo glm ‘m1’ y para el modelo bioclim ‘bc’, para tres registros con valores para las variables: cld6190_ann, dtr6190_ann, frs6190_ann, h_dem, pre6190_l1, pre6190_l10, pre6190_l4, pre6190_l7, vap6190_ann
cld6190_ann = c(40, 60, 80)
dtr6190_ann = c(40, 100, 140)
frs6190_ann = c(50, 125, 200)
h_dem = c(1000, 3000, 6000)
pre6190_l1 = c(50, 100, 150)
pre6190_l10 = c(50, 100, 150)
pre6190_l4 = c(50, 100, 150)
pre6190_l7 = c(50, 100, 150)
vap6190_ann = c(50, 100, 250)
pdt = data.frame ( cbind (cld6190_ann, dtr6190_ann, frs6190_ann, h_dem, pre6190_l1, pre6190_l10, pre6190_l4, pre6190_l7, vap6190_ann))
pdtpredict (m1 , pdt) #categorias de prediccion= 1: menos frecuente encontrar especies, 3: mas frecuente encontrar especies## 1 2 3
## -0.008807341 0.525429699 0.984308281predict (bc , pdt)7.3.3.1 Graficos de respuesta para cada variable
response (bc)
7.3.3.2 Mapa de puntajes de idoneidad
En el espectro de colores del mapa se observan los valores mas altos como aquellos donde es idoneo encontrar estas especies, al contrario en el espectro, donde es menos idoneo
p <- predict (predictors,m1)
plot (p)